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Background 

The  purpose  of  this  paper  is  to  report  the  results  of  testing  the  Fast 
Hartley  Transform  (FHT)  and  comparing  it  with  the  Fast  Fourier 
Transform  (FFT).  All  the  definitions  and  equations  in  this  paper  are 
quoted  and  cited  from  the  series  of  references.  The  author  of  this  report 
developed  a  Fortran  program  which  computes  the  Hartley  transform. 
He  tested  the  program  with  a  generalized  electromagnetic  pulse 
waveform  and  verified  the  result  with  the  known  value. 


1.  Introduction 

Fourier  analysis  is  an  essential  tool  to  obtain  frequency  domain 
information  from  transient  time  domain  signals.  The  FFT  is  a  popular 
tool  to  process  many  of  today's  audio  and  electromagnetic  signals. 
System  frequency  response,  digital  filtering  of  signals,  and  signal 
power  spectrum  are  the  most  practical  applications  of  the  FFT. 
However,  the  Fourier  integral  transform  or  the  FFT  requires  the 
computer  resources  appropriate  to  the  complex  arithmetic  opera¬ 
tions.  On  the  other  hand,  the  FHT  can  accomplish  the  same  results 
faster  and  requires  fewer  computer  resou rces^fTJpThe  FHT  is  twice  as 
fast  as  the  FFT,  uses  only  half  the  computer  resources,  and  so  could  be 
more  useful  than  the  FFT  in  typical  applications  such  as  spectral 
analysis,  signal  processing,  and  convolution  iikThis  paper  presents 
a  Fortran  computer  program  for  the  FHT  algorithm  along  with  a  brief 
description  and  compares  the  results  and  performance  of  the  FHT  and 
the  FFT  algorithms. 


2.  General  Description  of  FHT 


Equation  (1 )  defines  the  analytic  form  of  the  Hartley  transform,  and  (2) 
shows  its  inverse  transform,  which  switches  the  frequency  function 
back  into  the  time  domain  [1]: 


X(t)cas(2n ft)dt  , 


0) 


X(t)  = 


H(f)  cas  (2;i/f)#  . 


where  cas  (2  n ft)  =  cos  (2  nft)  +  sin  (2  nft). 


(2) 
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The  cas  function  was  introduced  by  R.  V.  L.  Hartley,  who  first 
proposed  the  Hartley  transform  in  1942  [3]. 


The  above  equations  are  very  similar  to  the  Fourier  transform,  equa¬ 
tion  (3),  and  its  inverse,  equation  (4)  [1]. 


Xiiyz-sW  dt  , 


(3) 


X(t)  =  J  F(f)eiWdf  ,  (4) 

where  e'2Kft  =  cos  (2  nft)  +  j  sin  (2  nft )  and  e'l2%ft  =  cos  (2  nft)  -  j  sin  (2  nft). 
(These  are  known  as  Euler's  formulas.) 

Note  the  electrical  engineering  convention  of  labeling  the  imaginary 
unit  i  as  j  ( i  equals  the  square  root  of  -1). 


These  four  equations  deal  only  with  continuous  time  variables.  In  the 
real  world,  however,  signals  are  sampled  at  discrete  intervals  of  time. 
So  there  are  discrete  transforms  that  approximate  the  Fourier  integral 
and  the  Hartley  integral.  The  discrete  forms  of  the  Hartley  transform 
(DHT)  pairs  are  [1] 

N- 1 

H(f)  cas  (2  nft/N)  ,  (5) 

N  o 

and  V-, 

XU)  =  X  H<fi  cas  <2  xfl/N)  .  (6) 

f=  0 


The  discrete  Fourier  transform  (DFT)  pairs  are  very  similar  to  the  DHT 

[1J. 

N-l 

F(f)  =  lJJJ  X(t)e-i^f‘lN  ,  (7) 

,  N  (  =  0 

and  N-l 

XU)  =  £  F(/)c<W>W  ,  (8) 

f=  0 


As  can  be  seen  in  equation  (7),  the  DFT  requires  (N-l)  complex 
multiplications  and  (N-l)  complex  additions  to  compute  each  output 
point  (the  first  term  in  the  sum  involves  exp  (j  *  0)  =  1  and  therefore 
does  not  require  a  multiplication).  Thus,  to  compute  N  output  points, 
N(N  -  1)  complex  multiplications  and  the  same  number  of  complex 
additions  are  required.  Now  each  complex  multiplication  requires 
four  real  multiplications  and  two  real  additions.  Hence,  for  computa- 
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tion  of  all  the  output  points  for  N-point  data,  the  DFT  requires  4 N(N 
- 1)  real  multiplications  and  4N(N  - 1)  real  additions.  To  overcome  this 
computational  requirement,  the  FFT  algorithm  (Cooley  and  Tukey 
algorithm)  was  developed  for  the  machine  computation  of  a  complex 
Fourier  transform  [4].  The  FFT  uses  a  permutation  process  to  bisect  the 
data  sequence  until  data  pairs  are  reached.  The  fundamental  concept 
for  the  permutation  process  is  that  it  is  faster  to  divide  the  data  set  into 
pairs,  compute  the  transform  of  the  pairs  individually,  and  recombine 
them  to  make  the  entire  transform  rather  than  to  compute  the  trans¬ 
form  as  a  whole  data  set.  The  Fourier  transform  of  the  time  domain 
data  set  can  be  obtained  by  superimposing  all  permutated  data  pairs. 
An  N- point  FFT,  where  N  is  a  power  of  2,  requires  2N  log2  N  real 
multiplications  and  3 N  log2  N  real  additions,  which  means  a  factor  of 
about  200  times  less  multiplication  for  N  equals  1024  time  data  points 
[5].  Similarly,  Bracewell  developed  a  fast  algorithm  for  the  Hartley 
transform.  However,  to  use  the  algorithm,  the  decomposition  formula 
that  expresses  a  complete  DHT  in  terms  of  its  half-length  subse¬ 
quences  is  required.  Bracewell  has  shown  the  following  decomposi¬ 
tion  formula  by  application  of  the  shift  and  the  similarity  theorems  for 
the  DHT  [6].  Using  an  equivalent  concept  for  the  Fourier  transform, 
the  similar  decomposition  formula  can  be  defined  [1]. 

H(J)  =  Hx(f)  +  H2(f)  cos  (2  nf/Ns) 

+  (Ns  ~f)  sin  (2  nf/Ns)  ,  w 

F(f>  =  Fl(f)  +  F2(j)eJW"’  ,  (10) 

where  N  is  the  number  of  elements  in  the  half-length  sequence,  and 
thus  Ns  =  N/2  for  a  data  set  of  N  elements 

As  can  be  seen,  there  is  one  important  difference  between  equations 
(9)  and  (10).  While  the  FFT  decomposition  formula  is  symmetric,  the 
FHT  decomposition  formula  is  asymmetric  because  of  the  sine  coeffi¬ 
cients  on  both  H(/)  and  H(N  -f).  This  aysmmetric  matrix  processing 
requires  special  handling  such  as  retrograde  indexing  for  computer 
implementation.  The  retrograde  indexing  behavior  can  be  described 
by  "using  an  independent  variable  as  an  index  for  the  elements 
multiplied  by  the  sine  coefficients.  This  index  decreases  while  the 
other  indexes  increase  [1]." 


3.  Comparison 

The  major  difference  between  the  two  transform  algorithms  is  the  real 
function  cas  in  the  Hartley  transform  and  the  complex  exponential 
term  in  the  Fourier  transform.  Since  real  arithmetic  is  much  simpler 


than  complex  computation,  the  FHT  is  faster  than  the  FFT  and  requires 
fewer  computer  resources.  Furthermore,  the  complex  Fourier  spec¬ 
trum  can  be  obtained  from  the  Hartley  transform.  It  is  faster  to 
generate  the  Fourier  transform  and  power  spectrum  with  the  FHT 
than  with  the  FFT,  because  the  FHT  algorithm  uses  real  rather  than 
complex  quantities  and  so  requires  fewer  floating-point  operations. 

If  a  function  can  be  expressed  uniquely  into  even  and  odd  parts,  then 
from  the  even  and  odd  parts  the  original  function  can  be  uniquely 
reconstructed.  Based  upon  the  symmetrical  property  for  the  even 
function  and  the  asymmetrical  property  for  the  odd  function,  the 
following  relationships  can  be  established  [6]. 

Let  H(f)  =  £(/)  +  0(/),  where  E(f)  and  0(f)  are  the  even  and  odd  parts  of 
H(f),  respectively.  Then 

E(f)  =  =  I  v(t )  cos  2  nftdt  (11) 

and  2  J- 

0(f)  I  Vd)  sin  2  nftdt  .  (12) 

According  to  equations  (11)  and  (12),  the  Fourier  transform  can  be 
obtained  from  the  Hartley  transform  simply  by  reflections  and  addi¬ 
tions  of  the  even  and  odd  parts.  Since  the  real  part  of  the  FFT  is  equal 
to  the  even  part  of  the  FHT,  and  the  imaginary  part  of  the  FFT  equals 
the  negative  odd  part  of  the  FHT,  the  real  and  imaginary  parts  of  the 
FFT  can  be  obtained  from  the  FHT  according  to  the  following  equa¬ 
tions  [1]: 

Fr  =  H(f)  +  H(N  -f)  (13) 

Fim  =  Hif)  -  H(N  -f)  ,  (14) 

where  Fr  is  the  real  portion  of  the  complex  Fourier  transform,  Fm  is  the 
imaginery  portion,  and  N  is  the  number  of  elements  in  the  data  set. 

The  power  spectrum  can  be  obtained  directly  from  the  FHT  using  the 
following  equation  [1]: 

Ps(f)  =  [W2  +  H(N-tf\/ 2.  (15) 

where  Ps  is  the  power  spectrum. 

While  the  FFT  computes  the  square  of  the  real  and  imaginary  parts 
and  sums  the  two  values  at  a  given  frequency,  the  FHT  squares  and 
sums  the  two  values  of  the  Hartley  transform  at  the  positive  and 
negative  frequencies.  Since  the  FHT  computes  the  energy  content  in 
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the  positive  and  negative  frequency  domain,  a  factor  of  two  is  required 
in  the  above  equation. 

The  convolution  theorem  is  almost  identical  between  the  Hartley  and 
the  Fourier  transforms.  Equation  (16)  summarizes  the  convolution 
theorem  for  the  Hartley  transform,  corresponding  to  equation  (17)  for 
the  Fourier  transform  [1]. 

flit)  ©/2(f)  =  Hx{f)  HM  +  ,  (16) 

where  HJif)  is  the  even  part  of  Hff)  and  H^if)  is  the  odd  part  of  Hff). 

fiit)®m=Fi(f)Fif)  (17) 

The  ®  symbol  denotes  the  convolution  operation. 

Note  that  if  one  of  the  functions  being  convolved  is  either  even  or  odd, 
then  the  convolution  theorem  for  the  Hartley  transform  reduces  to  the 
particularly  simple  form  indicated  below  [1]: 

Mt)®Ut)  =  Hx(f)H2(j)  .  (18) 


4.  Performance 

The  Fortran  program,  FHT.FOR,  is  presented  in  appendix  A  as  devel¬ 
oped  from  the  basic  program  presented  by  Brace  well  [2].  Using  the 
generalized  double  exponential  waveform  as  a  typical  electromag¬ 
netic  pulse  (EMP)  electric  field,  a  time  domain  data  file  is  generated. 
The  generalized  double  exponential  EMP  electric  field  time  behavior 
is  given  by 

E(t)  =  5.25*1 04[exp( -4*1 06/)-  exp  (-4.76*1 08/)]  ,  (19) 

in  volts  per  meter,  where  t  is  in  seconds.  This  pulse  has  a  peak  value 
of  50  kV/m,  a  10-  to  90-percent  risetime  of  about  5  ns,  and  a  time  to 
half- value  of  about  200  ns  [7].  Figure  1  plots  the  time  waveform  of  this 
constructed  high-altitude  EMP  electric  field.  The  time  domain  wave¬ 
form  has  1024  equally  spaced  points  with  a  time  increment  of  1  ns.  The 
user  can  choose  any  number  of  data  points  consistent  with  the 
available  computer  memory  size,  but  must  be  a  power  of  2  for 
algorithm  simplicity  and  additional  execution  speed.  The  sample  data 
have  to  be  equally  spaced  for  the  FHT  algorithm.  If  more  than  1024 
data  points  are  required,  the  FHT.FOR  program  has  to  be  modified  to 
expand  the  size  of  the  arrays.  The  time  increment  can  be  controlled 
also,  and  1  ns  is  used  as  an  example.  As  the  FHT.FOR  program 
executes,  the  user  has  an  option  to  choose  a  format  for  the  transformed 
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output  data.  One  format  has  frequency  and  magnitude  only  while  the 
other  one  includes  frequency  and  the  real  and  imaginary  parts  of  the 
FFT  spectrum.  Figure  2  shows  the  FHT.FOR  output,  and  figure  3 
shows  the  FFT  output  as  generated  with  the  HOBOII  signal-process¬ 
ing  software  package  [8].  As  can  be  seen  from  figures  2  and  3,  the  FFT 
and  FHT  results  appear  to  be  identical.  If  drawn  on  one  graph,  they 
would  overlay.  The  real  and  imaginary  parts  of  the  frequency  spec¬ 
trum  obtained  from  the  FHT  algorithm  were  verified  with  the  FFT 
output  from  HOBOII.  An  inverse  FFT  applied  to  the  FHT  results 
reproduces  the  input  waveform  which  verifies  the  accuracy  of  the 
FHT  algorithm. 

As  can  be  noticed  in  figures  2  and  3,  both  transform  algorithms 
produce  a  nonphysical  behavior  at  higher  frequencies.  To  explore  this 
deviation  and  compare  these  transforms  with  the  original  double 
exponential  EMP  spectrum,  another  Fortran  program,  BODE.FOR, 
was  developed.  This  program  is  included  in  appendix  B.  The 
BODE.FOR  program  computes  the  double  exponential  frequency 
spectrum  using  the  analytic  form  of  Fourier  transform  for  the  differ¬ 
ences  of  two  exponentials  [7].  The  calculated  spectrum  is  shown  in 
figure  4,  according  to  the  output  of  BODE.FOR.  In  figure  5,  the  FHT 
and  FFT  are  superimposed  on  the  analytic  spectrum.  Both  transforms 
trace  the  calculated  spectrum,  except  for  discrepancies  at  the  high 


Figure  1.  Generalized 
EMP  waveform. 


Time  (s) 
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frequency  as  previously  mentioned.  Since  the  double  exponential 
pulse  has  a  risetime  of  about  5  ns,  there  is  an  early  time  slope 
discontinuity  between  the  origin  and  the  first  sample  data  point  for  a 
1-ns  sample  time  increment.  This  situation  can  be  improved  if  we  take 
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Figure  4.  Calculated 
double  exponential 
EMP  spectrum. 


Solid  line,  calculated  double  exponential  EMP  spectrum. 
Dotted  line,  FFT  and  FHT. 
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more  data  points  in  the  early  time,  which  requires  a  smaller  time 
increment.  Another  difference  between  the  analytic  frequency  spec¬ 
trum  and  that  obtained  with  the  FHT  and  FFT  transforms  is  that  the 
low  frequency  limit  is  undetermined.  In  order  to  obtain  the  low 
frequency  value  of  the  transformed  data,  the  time  domain  data  set  has 
to  have  a  much  larger  time  window,  which  would  require  more  time 
to  process  the  transform.  The  time  increment  could  also  be  increased, 
but  would  lead  to  a  '  <s  of  the  high-frequency  content.  According  to 
the  Nyquist  sampling  theorem,  200,000  time  domain  data  points  are 
required  with  a  time  increment  of  0.5  ns  to  generate  a  double  exponen¬ 
tial  frequency  spectrum  with  the  FHT  or  FFT,  which  would  be  similar 
to  figure  4.  The  Nyquist  sampling  theorem  requires  that  the  sampling 
rate  must  be  at  least  twice  the  highest  frequency  of  interest  in  the 
waveform  being  sampled  [9]. 

When  the  speed  of  two  transforms  is  compared,  theoretically,  the  FHT 
is  faster  than  the  FFT;  however,  it  is  hard  to  compare  the  run  time  of 
these  two  programs  for  the  transform  process  itself  because  the 
programs  used  in  this  effort,  FHT.FOR  and  HOBOII,  are  implemented 
differently.  Also,  for  all  practical  purposes  both  codes  process  the 
transform  within  a  few  seconds. 


5.  Conclusions 

Since  the  FHT  uses  only  real  valued  functions,  there  is  no  need  for 
complex  calculations,  which  implies  faster  run-times  and  less  com¬ 
puter  memory  to  process  a  signal  in  comparison  to  a  typical  FFT 
algorithm.  Finally,  the  FHT  uses  fewer  operations  to  transform  a  given 
signal,  so  there  are  fewer  round-off  errors. 

The  example  illustrated  does  not  prove  that  the  FHT  is  superior  to  the 
FFT,  although  it  demonstrates  that  the  FHT  is  fully  compatible  with 
the  FFT.  However,  if  large  amounts  of  data  are  being  manipulated, 
there  could  be  a  significant  difference  in  speed  and  the  amount  of 
memory  required.  The  FHT  offers  better  performance  us'  ig  fewer 
computational  resources. 

The  FHT.FOR  program  can  be  further  developed  to  perform  a  convo¬ 
lution  that  requires  only  real  arithmetic.  An  FHT  convolution  algo¬ 
rithm  would  be  faster  and  easier  to  implement  than  a  similar  algo¬ 
rithm  developed  for  the  FFT. 
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Appendix  A. — Fortran  Program  FHT.FOR 


APPENDIX  A 


$LARGE 

C  This  routine  generates  the  fast  Hartley  transform  (FHT) 

C  and  the  complex  form  of  fast  Fourier  transform  (FFD 
C  according  to  a  given  time  domain  function. 

REAL  L2 

INTEGER  P7/T/U,S2,Q,D,E,SS,SO/P,OPTION 

DIMENSION  F(10/1024)/R(1024),X(1024)/FR(1024)/RMAG(1024)  DIMENSION 
S(1024),C(1024)/M(11)/RR(1024)/XX(1024)  CHARACTER*20  FILN AM,FILEN AM 

C  Opening  the  output  files. 

WRITE(V)  TYPE  OUTPUT  FILENAME  FOR  THE  TIME  DOMAIN'  READ(V) 
FILNAM 

OPEN(3,FILE=FILNAM/STATUS='NEW') 

WRITEf*,*)  TYPE  OUTPUT  FILENAME  FOR  THE  FREQUENCY  SPECTRUM' 
READ(V)  FILENAM 

OPEN  (10,FILE=FILEN  AM,ST  ATUS='NEW') 

C  Initialization. 


WRITE(V)  'HOW  MANY  DATA  POINTS.#7 

WRITE(*,*)  TOWER  OF  2  BUT  LESS  OR  EQUAL  TO  1024'  READ(*,*)  NU 
P  =  0 

1  NUM  =  NU/2 
P  =  P  +  1 
NU  =  NUM 

IF  (NUM.  GE  .2)  GOTO  1 
N4  =  2**(P-2) 

N2  =  N4  +  N4 
N  =  N2  +  N2 
N7  =  N  - 1 
P7  =  P  - 1 

WRITE(  V)  'WHAT  IS  YOUR  TIME  INCREMENT  ?'  READ(V)  DTM 
WRITEf*,*)  'WHAT  ARE  YOUR  PARAMETERS  ?'  WRITE(V)  'E(t)  =  EO*(EXP(- 
AP*t)-EXP(-BT*t)'  WRITE(*,*)  'EO=?,  AP=?,  BT=?' 

READ(*,*)  EO,AP,BT 
TM  =  0.0 
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C  Generating  the  time  domain  waveform. 

DO  10 1  =  0,N7 

FNF  =  EO*(EXP(-APTM)-EXP(-BT*TM))  WRITE(3,*)  TM,FNF 
TM  =  TM  +  DTM 
DF  =  1.0/((N-1)*DTM) 

FR(I)  =  PDF 
F(0,D  =  FNF 
F0,I)  =  FNF 
10  CONTINUE 

C  Generating  the  power  of  2  numbers. 


1  =  1 

M(O)  =  1 
M(l)  =  2 

20  M(l+1)  =  M(I)  +  M(I) 

1  =  1  +  1 

IF  (I.  LT  .P)  GOTO  20 


C  Get  the  sin  and  cos  coefficients. 

PI  =  4*ATAN(1.0) 

W  =  2*PI/N 
A  =  0 

DO30I  =  l,N 
A  =  A  +  W 
S(I)  =  SIN(A) 

C(I)  =  COS(A) 

30  CONTINUE 

C  Start  permutation. 

J  =  -l 
I  =  -l 
50  1  =  1  +  1 
T  =  P 
40  T  =  T-l 

J  =  }  -  M(T) 

IF  (J.  GE  .-1)  GOTO  40 
J  =  J  +  M(T+1) 

IF(I.LE.p  GOTO  50 
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T  =  F(0,I+1) 

F(0,I+1)  =  F(OJ+l) 

F(OJ+l)  =  T 
IF  (I.  LT  .N-3)  GOTO  50 

C  First  stage. 

DO  60 1  =  0,N-2,2 

FO/I)  =  F(0,D  +  F(0,I+1)  F(U+1)  =  F(0,I)  -  F(0,I+1) 
60  CONTINUE 

C  Second  stage. 


U  =  P7 
SS  =  4 

DO  90  L  =  2,P7 
S2  =  SS  +  SS 
U  =  U-1 
SO  =  M(U-l) 

DO  100  Q  =  0,N7,S2 
I  =  Q 
D  =  I  +  SS 

F(L+1,I)  =  F(L,I)  +  F(L,D)  F(L+1,D)  =  F(L,I)  -  F(L,D)  K  =  D  - 1 
DO  110  J  =  SO^SO 
1  =  1+1 
D  =  I  +  SS 
E  =  K  +  SS 

Y  =  F(L,D)*C(J)  +  F(L,E)*S(J)  Z  =  F(L,D)*S(J)  -  F(L,E)*C(J)  F(L+I,I)  =  F(L,I) 
+  Y 

F(L+1,D)  =  F(L,I)  -  Y 
F(L+1,K)  =  F(L,I)  +  Z 
F(L+1,E)  =  F(L,K)  -  Z 

K  =  K  +  1 
110  CONTINUE 

E  =  K  +  SS 

100  CONTINUE 

SS  =  S2 

90  CONTINUE 

C  Normalizing  the  fast  Hartley  transform's  magnitude. 

RMAGfO)  =  F(L,0)/(DF*N) 
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RR(0)  =  F(L,0)/(DF*N) 

XX(O)  =  0 


C  Select  either  ma#nit#de  of  the  FHT 
C  or  the  real  and  #mag#nary  part  of  the  FFT. 

WRITE/,*)  'SELECT  THE  DESIRED  OUTPUT'  WRITE/,*)  T  =  MAGNITUDE,  2  = 
REAL  &  IMAGINARY'  READ/,*)  OPTION 
IF  (OPTION.  EQ  .1)  THEN 

C  The  fast  Hartley  transform's  magnitude  is  derived  from  the 
C  power  spectrum 

DO  120 1  =  N-l,(N/2)+l,-l 

F(L,I)  =  F(L,I)/(N*DF) 

F(L,N-I)  =  F(L,N-I)/(N*DF) 

RMAG(I)  =  SQRT((F(L,I)**2  +  F(L,N-I)**2)/2)  WRITE  (10/)  FR(N-I),  RMAG(I) 
120  CONTINUE 

C  Get  the  real  and  imaginary  parts  of  the  fast  Fourier 
C  transform  from  the  fast  Hartley  transform. 

ELSE 

DO  130  I  =  N-l,(N/2)+l,-l 
B  =  F(L,N-I) 

RR(I)  =  F(L,I)  +  B 
XX(I)  =  F(L,I)  -  B 
IF  (I.  LE  .N-I)  THEN 
J  =  I 

ELSE 

J  =  N  - 1 
RR(I)  =  RR(J) 

XX(J)  =  -XX(J) 

ENDIF 

R(I)  =  RR(I)/(2*N*DF) 

X(I)  =  XX(I)/(2*N*DF) 

WRITEdO/)  FR(N-I),R(I),X(I) 

130  CONTINUE 

ENDIF 

END 
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Appendix  B. — Fortran  Program  BODE.FOR 


APPENDIX  B 


C  This  is  a  program  to  compute  the  analytic  frequency  spectrum 
C  of  a  double  exponential  pulse. 

C  Use  COMPLEX  constants,  variables,  operations,  and  functions. 

C  Input  variables  (all  real) : 

C  First :  Starting  frequency 

C  Last  :  Stopping  frequency 

C  INC  :  Additive  frequency  increment 

C  Intermediate  variables,  all  complex 

C  D1 :  First  factor  in  denominator 
C  D2  :  Second  factor  in  denominator 

C  These  factors  are  set  up  for  using  the  CMPLX  function,  which 
C  converts  from  the  form  of  two  REAL  values,  representing  the 
C  real  and  imaginary  parts  of  the  complex  numbers,  to  the  form 
C  of  one  Fortran  COMPLEX  number. 


REAL*8  K,  FIRST,  LAST,  INC,  OMEGA,  ABSVAL,  FRE  COMPLEX*!  6  E,  D1,D2 
CHARACTER*20  FILNAM 

C  Opening  the  output  file. 

WRITE(*,*)  TYPE  OUTPUT  FILENAME' 

READ(*,*)  FILNAM 

OPEN(3,FILE=FILNAM,STATUS='NEW') 

C  Read  parameters,  validate. 

WRITE(*,*)  'FIRST  =  ?,  LAST  =  ?,  INC  =  ?'  READ(*,*)  FIRST,  LAST,  INC 
IF  (FIRST.  GE  .LAST)  THEN 

WRITE(V)  'INVALID  DATA;  PROGRAM  ABORTED'  STOP 
END  IF 

C  Set  frequency  to  starting  value. 

OMEGA  =  FIRST*8*ATAN(LO) 

C  While  OMEGA  <=  LAST 

10  CONTINUE 

IF  (OMEGA.  LE  .LAST)  THEN 
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D1  =  CMPLX(4.0D6,0MEGA) 

D2  =  CMPLX(4.76D8,0MEGA) 

E  =  2.47D13/(D1*D2) 

C  Get  complex  absolute  value,  =  magnitude  of  output. 
ABSBAL  =  CDABS(E) 

FRE  =  OMEGA/ (8*ATAN(1.0))  WRITE (3,*)  FRE,ABSVAL 
C  Incrementing  frequency. 


FRE  =  INC  +  FRE 

OMEGA  =  FRE*8*  AT  AN  (1.0) 

GOTO  10 

END  IF 

END 
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